/*
 * Class:        BrownianMotion
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
package umontreal.ssj.stochprocess;
import umontreal.ssj.rng.*;
import umontreal.ssj.probdist.*;
import umontreal.ssj.randvar.*;

/**
 * This class represents a *Brownian motion* process @f$\{X(t) : t \geq0
 * \}@f$, sampled at times @f$0 = t_0 < t_1 < \cdots< t_d@f$. This process
 * obeys the stochastic differential equation
 * @anchor REF_stochprocess_BrownianMotion_eq_Brownian_motion
 * @f[
 *   dX(t) = \mu dt + \sigma dB(t), \tag{Brownian-motion}
 * @f]
 * with initial condition @f$X(0)= x_0@f$, where @f$\mu@f$ and @f$\sigma@f$
 * are the drift and volatility parameters, and @f$\{B(t),  t\ge0\}@f$ is a
 * standard Brownian motion (with drift 0 and volatility 1). This process has
 * stationary and independent increments over disjoint time intervals (it is
 * a Lévy process) and the increment over an interval of length @f$t@f$ is
 * normally distributed with mean @f$\mu t@f$ and variance @f$\sigma^2 t@f$.
 *
 * In this class, this process is generated using the sequential (or random
 * walk) technique: @f$X(0)=x_0@f$ and
 * @anchor REF_stochprocess_BrownianMotion_eq_Brownian_motion_sequential
 * @f[
 *   X(t_j) - X(t_{j-1}) = \mu(t_j - t_{j-1}) + \sigma\sqrt{t_j - t_{j-1}} Z_j \tag{Brownian-motion-sequential}
 * @f]
 * where @f$Z_j \sim N(0,1)@f$.
 *
 * <div class="SSJ-bigskip"></div><div class="SSJ-bigskip"></div>
 */
public class BrownianMotion extends StochasticProcess {
    protected NormalGen    gen;
    protected double       mu,
                           sigma;
    // Precomputed values for standard BM
    protected double[]     mudt,
                           sigmasqrdt;

   /**
    * Constructs a new `BrownianMotion` with parameters @f$\mu=@f$ `mu`,
    * @f$\sigma=@f$ `sigma` and initial value @f$X(t_0) =@f$ `x0`. The
    * normal variates @f$Z_j@f$ in (
    * {@link REF_stochprocess_BrownianMotion_eq_Brownian_motion_sequential
    * Brownian-motion-sequential} ) will be generated by inversion using
    * `stream`.
    */
   public BrownianMotion (double x0, double mu, double sigma,
                          RandomStream stream) {
        this (x0, mu, sigma, new NormalGen (stream));
    }

   /**
    * Constructs a new `BrownianMotion` with parameters @f$\mu=@f$ `mu`,
    * @f$\sigma=@f$ `sigma` and initial value @f$X(t_0) =@f$ `x0`. Here,
    * the normal variate generator  @ref umontreal.ssj.randvar.NormalGen
    * is specified directly instead of specifying the stream and using
    * inversion. The normal generator `gen` can use another method than
    * inversion.
    */
   public BrownianMotion (double x0, double mu, double sigma, NormalGen gen) {
        this.mu    = mu;
        this.sigma = sigma;
        this.x0    = x0;
        this.gen   = gen;
    }


   public double nextObservation() {
        double x = path[observationIndex];
        x += mudt[observationIndex]
             + sigmasqrdt[observationIndex] * gen.nextDouble();
        observationIndex++;
        path[observationIndex] = x;
        return x;
    }

/**
 * Generates and returns the next observation at time @f$t_{j+1} =@f$
 * `nextTime`. It uses the previous observation time @f$t_j@f$ defined
 * earlier (either by this method or by <tt>setObservationTimes</tt>), as
 * well as the value of the previous observation @f$X(t_j)@f$. *Warning*:
 * This method will reset the observations time @f$t_{j+1}@f$ for this
 * process to `nextTime`. The user must make sure that the @f$t_{j+1}@f$
 * supplied is @f$\geq t_j@f$.
 */
public double nextObservation (double nextTime) {
        // This method is useful for generating variance gamma processes
        double x = path[observationIndex];
        double previousTime = t[observationIndex];
        observationIndex++;
        t[observationIndex] = nextTime;
        double dt = nextTime - previousTime;
        x += mu * dt + sigma * Math.sqrt (dt) * gen.nextDouble();
        path[observationIndex] = x;
        return x;
    }

   /**
    * Generates an observation of the process in `dt` time units, assuming
    * that the process has value @f$x@f$ at the current time. Uses the
    * process parameters specified in the constructor. Note that this
    * method does not affect the sample path of the process stored
    * internally (if any).
    */
   public double nextObservation (double x, double dt) {
        x += mu * dt + sigma * Math.sqrt (dt) * gen.nextDouble();
        return x;
    }


   public double[] generatePath() {
        double x = x0;
        for (int j = 0; j < d; j++) {
            x += mudt[j] + sigmasqrdt[j] * gen.nextDouble();
            path[j + 1] = x;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }

/**
 * Same as generatePath(), but a vector of uniform random numbers must be
 * provided to the method. These uniform random numbers are used to generate
 * the path.
 */
public double[] generatePath (double[] uniform01) {
        double x = x0;
        for (int j = 0; j < d; j++) {
            x += mudt[j] + sigmasqrdt[j] * NormalDist.inverseF01(uniform01[j]);
            path[j + 1] = x;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }


   public double[] generatePath (RandomStream stream) {
        gen.setStream (stream);
        return generatePath();
    }

/**
 * Resets the parameters @f$X(t_0) = \mathtt{x0}@f$, @f$\mu= \mathtt{mu}@f$
 * and @f$\sigma= \mathtt{sigma}@f$ of the process. *Warning*: This method
 * will recompute some quantities stored internally, which may be slow if
 * called too frequently.
 */
public void setParams (double x0, double mu, double sigma) {
        this.x0    = x0;
        this.mu    = mu;
        if (sigma <= 0)
           throw new IllegalArgumentException ("sigma <= 0");
        this.sigma = sigma;
        if (observationTimesSet) init(); // Otherwise not needed.
    }

   /**
    * Resets the random stream of the normal generator to `stream`.
    */
   public void setStream (RandomStream stream) { gen.setStream (stream); }

   /**
    * Returns the random stream of the normal generator.
    */
   public RandomStream getStream() { return gen.getStream (); }

   /**
    * Returns the value of @f$\mu@f$.
    */
   public double getMu() { return mu; }

   /**
    * Returns the value of @f$\sigma@f$.
    */
   public double getSigma() { return sigma; }

   /**
    * Returns the normal random variate generator used. The
    * @ref umontreal.ssj.rng.RandomStream used by that generator can be
    * changed via `getGen().setStream(stream)`, for example.
    */
   public NormalGen getGen() { return gen; }


    // This is called by setObservationTimes to precompute constants
    // in order to speed up the path generation.
   protected void init() {
        super.init();
        mudt       = new double[d];
        sigmasqrdt = new double[d];
        for (int j = 0; j < d; j++) {
            double dt     = t[j+1] - t[j];
            mudt[j]       = mu * dt;
            sigmasqrdt[j] = sigma * Math.sqrt (dt);
        }
     }

}